Critical velocity, vortex shedding and drag in a unitary Fermi superfluid 
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We study the real-time motion of a microscopic object in a cold Fermi gas at unitary conditions by 
using an extended Thomas- Fermi density functional approach. We find that spontaneous creation 
of singly quantized vortex-antivortex pairs occurs as a critical velocity is exceeded, which leads to a 
drag between the moving object and the Fermi gas. The resulting force is linear in the velocity for 
subsonic motion and becomes quadratic for supersonic motion. 
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I. INTRODUCTION 

A Fermi gas of atoms at unitary conditions is charac- 
terized by a divergent s-wave scattering length , which 
makes it a unique system, being at the same time di- 
lute and strongly interacting. In this regime all scales 
associated with interactions disappear from the problem 
and the energy of the system is expected to be propor- 
tional to that of a non interacting fermions system. The 
unitary Fermi gas (UFG) at low energies is known to 
be superfluid, and the most striking manifestation of its 
phase coherence has been the experimental observation 
of vortices in a rotating UFG of ^Li atoms • 

Vortex dipoles in a Fermi gas were studied theoreti- 
cally, within an extended Thomas-Fermi approach, in the 
BCS limit [3]. Vortex solutions within the same approach 
have been studied in Ref. (j], where the rotation of the 
Fermi system was predicted to generate a giant vortex 
in the presence of strong anharmonicity in the confining 
potential. Recently f5| a time-dependent Bogoliubov de- 
Gennes (BdG) approach has been used to study the 3-D, 
real-time formation of vortices in a UFG. Surprisingly, 
one of the main conclusion in Ref. [Sj is that the sys- 
tem remains superfluid even when stirred at supercritical 
speed. Single vortex solutions in the Ginzburg-Landau 
regime of a trapped superfluid Fermi gase have been stud- 
ied in Ref. and Ref. 0. In a dilute Fermionic super- 
fluid a vortex state is characterized by a strong density 
depletion along the vortex core. The depletion is however 
not complete, according to the BdG calculations of Refs. 
ii- 

Recently it has been remarked [l0| that the superfluid 
unitary Fermi gas can be efficiently described at zero tem- 
perature by phenomenological density functional theory. 
Density functionals of different flavours have been pro- 
posed by different theoretical groups. Bulgac and Yu 
have introduced a superfluid density functional based on 
a Bogoliubov-de Gennes approach to superfluid fermions 
[ill [l2| . Papenbrock and Bhattacharyya [l^ have in- 
stead proposed a Kohn-Sham density functional with an 
effective mass to take into account nonlocality effects. 
Here we adopt instead the extended Thomas-Fermi func- 
tional of the UFG that we have proposed few years ago 



[l3 | and which has been used recently to successfully ad- 
dress a number of properties of such system [T^ - fTsj . 

Our extended Thomas- Fermi (ETF) density functional 
approach for the UFG [l3, [ill is based on the following 
extended hydrodynamics equations: 
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where n(r, t) is the time-dependent scalar density field 
and v(r,i) the time-dependent velocity field. Here 
, t)_is the external potential and ^ = 0.40 and X — 1/4 
The above equations describe accurately vari- 
ous static and dynamical properties of the UFG. The 
term multiplied by A is found to be crucial to accurately 
describe surface effects, in particular in systems with a 
small number of atoms, where the Thomas-Fermi (local 
density) approximation fails [l3|- We have shown that 
when fast dynamical processes occur and/or when shock 
waves come into play such term is necessary also in the 
large N limit 18| , where results in quantitative agreement 



with experiments of real-time collision of strongly inter- 
acting Fermi gas clouds at unitarity fl^ can be obtained 

We use this method here to study the motion of a mi- 
croscopic 2-dimensional object with circular shape in a 
UFG and the associated process of vortex shedding, once 
the object velocity exceeds a critical value. Experimen- 
tal realization of this geometry would imply, for instance, 
moving a far-detuned laser beam through a trapped con- 
densate. 



II. CRITICAL VELOCITY AND DRAG IN THE 
UNITARY FERMI GAS 

For a dilute Fermi gas the long wavelength elementary 
excitations are sound waves, and the Landau criterion 
for the critical velocity Vc = {t/p)min for breakdown of 
superfluidity gives Vc = ci. However, for fluid flowing 
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past an object (or an atomic impurity as well), the local 
velocity near the object surface can become supersonic 
even when the fluid velocity far from the object, foo, is 
subsonic (the maximum velocity of, e.g., a 2-dimensional 
flow of an incompressible fluid past a circular obstacle 
is reached on the perimeter of the obstacle, where it is 

An estimate for the critical velocity for a UFG can be 
obtained as follows. At stationary conditions Eq. ([2]) 
provides the Bernoulli equation for the UFG: (in the fol- 
lowing a ^^£(37r2)2/3) 

-X^^^l^ + an^/' + Uir) + ^v' = const (3) 
2m y/n 2 

By assuming that the quantum term in the previous 
equation is negligible in the spirit of the long-wavelength 
limit approximation, and calling no the (uniform) density 
far from an impenetrable object, where v ~ Woo, one 
flnds: 

n{r)^[nl^' + ^{vl-v')r/' (4) 

outside the object and n = within it (the interaction 
U between the object and the fluid only provides the 
excluded volume condition). 

The above equation may be written in terms of the 
local 

cioc = \I^VF{r) = ^2an{vY/^/{im) (5) 

and bulk 

Cl = \j^VF,oo = \/2anl^^/{3m) (6) 

speeds of sound as: 

cL = cl + (vl - v^)/3 (7) 



By combining the two equations ([TJ and ([2]) one can de- 
rive the equation for the momentum conservation (sum- 
mation over repeated indices is implied): 

dtJk + d^T.k + pdk [U/m) = (9) 

where Ju = pvk is the supercurrent density and the stress 
tensor Ti^ is defined as: 

T,, ^ pv^Vk + <5,fc(^p5/3) ~ h^fpd^dkHp) (10) 
5m 4 2m 

(to derive the above equation the following identity has 
been used: p~^di{pdidkln{p)) = 2^k{^^^^y/p/ ^/p)). 

Note that although in a superfluid there is no frictional 
viscosity, nonetheless shear stress may arise from density 
(pressure) gradients (third term in Eq. (|10p . As a conse- 
quence, vortex formation and drag are possible even with 
no viscosity. 

The force exerted on a condensate by an obstacle mov- 
ing through it can be calculated from the rate of momen- 
tum transfer to the fluid. By integration, one finds the 
drag force (per unit mass): 

Fk^dt [ dQJk = - / dSn,T,fc - / dilpdhU (11) 
Jo is Jn 

Here S is the object external surface and n is a unit 
vector directed along the outward normal. In the case of 
an homogeneous fiow past an impenetrable object only 
the first term contributes. In the present case, however, 
where a partially penetrable object is used (see the fol- 
lowing) both contributions are present. 

In our calculations we used the full expression (|TOl) for 
the stress tensor to compute the drag ((TT|) . Notice that 
the quantum term in (jlOp is expected to be negligible 
only when Vao ^ {h/mR), where R is the diameter of 
the moving object: this is not the case here since these 
two terms are of comparable magnitude. 

III. METHODS AND CALCULATIONS 



When V ^ Qoc local instabilities develop, leading to the 
release of vortices [l^l- The maximum local velocity of 
a fluid flowing past an impenetrable cylindrical object, 
normally to its axis, is u ^ 2voo (on the surface of the 
cylinder and tangent to it). Using this value in ([7]) one 
can solve the equation for Uqo , thus providing an approx- 
imate value for the critical velocity of the fluid flowing 
past a stationary object (or, equivalently, of a moving 
object in the fluid at rest): 

vc = ci/Vb (8) 



This value is very similar to the one, Vc ~ ^/2/l^Cl 
obtained for the BEG case, using a similar argument, in 
Ref. 20]. For a spherical object the equatorial velocity 
is instead v ~ 3woo/2, giving for the critical velocity of a 
sphere moving in a UFG Vc = \/3/8ci. 



By using a Madelung transformation, equations ([T]) 
and ([2) can equivalently be written in the form of a 
time-dependent nonlinear Schrodinger equation (NLSE) 
[l3 | involving the complex order parameter ^{r,t) ~ 

d 

i/i— * = H"^ (12) 
dt ^ ' 

where 

7?= -^V2 + 2C/(r- Vt) + 2a|*|^/3 (13) 
4m 

The link between the two descriptions is provided by 
the deflnition 

v(r,i) = ^V0(r,i) (14) 
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for the velocity field associated with the phase of the 
order parameter 9{r,t). V in Eq. is the (constant) 
velocity of the moving object. 

By means of a Galilean transformation to the reference 
frame moving with the object (which will thus appear as 
stationary in our simulations), the NLSE to be solved 
becomes: 



H + iKV-V 



(15) 



For simplicity, we model the microscopic object inside 
the moving Fermi gas by means of a repulsive cylindri- 
cal "wall" and consider the flow motion perpendicular 
to it. Due to translational invariance along the axis of 
the cylinder (which we define as the z axis) , the problem 
thus reduces to finding the density and the velocity of 
the fluid in the (x, y) plane. 

We have numerically solved this 2-D NLSE equation to 
obtain the long-time dynamics of the fluid due to the mo- 
tion of the object moving along the x direction. From the 
time-dependent calculated density n(x,y,t) and velocity 
v(a:, J/, t) we computed the drag acting on the object, ac- 
cording to pT|) . 

Our simulation region is a square cell of side a = 
3.2 fim. We used a uniform mesh to represent the wave- 
function ^ff, on 512 X 512 points in the {x,y) plane. We 
have used the Runge-Kutta-Gill fourth-order method [2l| 
to propagate in time the solutions of the NLSE. To accu- 
rately compute the spatial derivatives appearing in our 
NLSE, we used a 13-point flnite-difference formula p^ . 
To avoid the outgoing sound waves and/or emitted vor- 
tices to interfere with the fluid dynamics after being re- 
flected on on the grid boundaries, we use an exponential 
absorbing buffer located in the periphery of the cell, as 
described in Ref . [l^ . The waves can travel freely in the 
undamped region (which occupy most of the simulation 
cell) but are quickly attenuated as they enter the exter- 
nal region, thus preventing unwanted interferences which 
might spoil the results. 

In the following we take m equal to the mass of a ^Li 
atom. We consider two different density values ("high 
density" and "low density" in the following) for the UFG, 
namely uq — 8600 /im^'^ and hq = 880 /im~^, corre- 
sponding to interparticle distances d ^ 5 x 10~^ /im and 
d ^ 0.1/im, respectively. 

At equilibrium, the repulsive cylindrical "wall" results 
in the formation of a circular cavity void of atoms. The 
object potential U and the associated initial density pro- 
flle are shown in Fig. ([1]). 



IV. RESULTS AND DISCUSSION 

We made a series of time-dependent calculations, solv- 
ing Eq. (|15p for different values of the speed V, and 
following the dynamics of the systems for several /xs. We 
flnd that there is a critical value Vc (which will be quanti- 
fled in the following) for the object speed V which sepa- 
rates two distinct regimes. When V < Vc the fluid profile 
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FIG. 1: Equilibrium density profile at t=0. The dotted line 
represents (in arbitrary units) the repulsive potential due to 
the object. 



rapidly evolves with time into a stationary configuration. 
Both the density and the velocity field for such final con- 
figuration are fore-to-aft symmetric. This implies that 
the drag force on the object is zero, which is another ver- 
sion of the well-known D'Alembert paradox in classical 
fluids. 

In Fig. ([2]) (lower curve) we show the calculated drag 
force Fx (computed using Eq. (ITT|) ') as a function of time 
when V < Vc- After a transient where the fluid acco- 
modates under the sudden acceleration of the object, the 
drag goes eventually to zero when the symmetric (sta- 
tionary) conflguration is reached. 

Above Vc, however, the Fermi gas dynamics changes 
dramatically: linear (antilinear) vortices are sponta- 
neously created in pairs close to the surface of the object, 
together with the emission of sound waves. 

We show in Fig. ^ a sequence of images taken dur- 
ing the real-time evolution of the system when V > Vc- 
The object is moving from left to right. First a localized 
bow wave moving with supersonic velocity is emitted in 
front of the cylinder and rapidly moves ahead. This is 
the result of a " shock" wave produced by the sudden ac- 
celeration of the object (we recall that the initial state 
is with fluid at rest and the object moving at constant 
velocity). 

Then a vortex pair is emitted in the rear, trailing be- 
hind (meaning that the pair velocity is lower than that of 
the moving object). The initial pair is soon overtaken by 
a pair emitted successively: the vortex lines apparently 
form a temporarily bound state (but they will eventually 
split at later times, not shown). 

As shown in the upper curve of Fig. ([2]), in this case 
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FIG. 2: Drag force during the real-time evolution of the sys- 
tem. Upper curve: V > Vc, lower curve: V < Vc- 




FIG. 3: From left to right, top to bottom: configurations 
at increasing times for V/ci = 0.76. Only a portion of the 
simulation cell is displayed. 

the instantaneous calculated drag force relaxes, after a 
transient, towards a finite nonzero value (with oscilla- 
tions that reflect the quasi-periodic emission of vortex 
pairs) . 

In the simulations of Fig. (jSj the velocity V is greater 
than Vc, but still below the speed of sound ci. 

The sequence shown in Fig. Q is instead obtained for 




FIG. 4: From left to right, top to bottom: configurations at 
increasing times for velocity V/ci = 1.44. Only a portion of 
the simulation cell is displayed. 

a supersonic motion of the object, y/ci = 1.44. It ap- 
pears that the vortex shedding frequency is considerably 
increased with respect to the previous case (similarly to 
what happens for a BEC, where the shedding frequency is 
[2^ 1 oc y^. Vortex- antivortex pairs are emitted in a semi- 
continuous way on different part of the rear section of the 
object, leading to parallel rows of connected vortical lines 
that eventually decay into separated vortices. Note also 
the appearance of a rather structured bow sound wave 
pattern moving along with the object. As time proceeds, 
the wake region behind the object becomes turbulent due 
to the superposition of more and more vortices and sound 
waves. 

In order to better characterize the vortex structure, we 
have selected (from another simulation where V is just 
above Vc-, so that vortex pairs are emitted with low fre- 
quency and thus are well separated from one another) 
a configuration after a pair of vortex lines have been 
emitted and moved away from the object, and closely 
analyzed the vortex structure. We find that the vor- 
tex is singly quantized, and that the two vortices of the 
pair have opposite polarization. The velocity field (not 
shown) follows very closely the ideal vortex velocity pro- 
file, v(r) = h/(2mr) (here r = a/x^ 4- y^). 

The vortex line structure has an empty core, as shown 
in Fig. ([5]), while the core size scale is set by kp^. This 
is in contrast with calculations based on Bogoliubov-De 
Gennes calculations [1, Q where a partially filled core 
(between 0.2 and 0.3 no) is predicted due to the presence 
of some normal liquid coming from pair breaking in the 
core region, where the velocity is higher. We will discuss 
in the following the reason for such discrepancy. 
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FIG. 5: Density profile in the vortex core region. 
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FIG. 6: Current density distribution near the vortex core. 
Squares: our results. Solid line: BdG calculations from Ref. 



We computed the superfluid current density j — p'v cir- 
culating around the vortex core (see Fig. (jG])), and found 
a peak value at rmax, in a remarkable overall agreement 
with the T=0 Bogoliubov-De Gennes calculations of Ref. 
[§| (shown in Fig. ([6]) with a solid line). 

We find the agreement with the results of Ref. 9] par- 
ticularly rewarding since it shows that an important su- 
perfluid observable related to vorticity can indeed be de- 



scribed accurately by our Density Functional approach. 

The scale of the maximum circulating current is set by 
the critical velocity which is determined by pair-breaking 
on the BCS side and by collective excitations on the BEG 
side 0. This is why our approach, which cannot obvi- 
ously describe single-particle processes like pair-breaking, 
is nonetheless able to reproduce the current pattern. For 
T < Tmax the kinetic energy cost associated with the cur- 
rent flow becomes larger than the condensation energy. 
Thus Tmax givcs an estimate of the distance from the core 
center below which superfluidity is partially suppressed. 

We note that the peak position r,nax in Fig- & coin- 
cides with the value of r where the vortex velocity field 
becomes equal to the local sound velocity, i.e. when 



j2an(r)^/^/3m (16) 

2mr V 

where n{r) is the vortex core density profile shown in Fig. 

In Ref. [16j we have found that, at unitarity, our ap- 
proach leads to a maximum Josephson current across a 
barrier which is practically the same as the one obtained 
from BdG. This suggests that, at unitarity, the current 
is limited by Landau's criterion for the creation of col- 
lective excitations, and not by single particle excitations. 
It is then not surprising that also in the present case our 
maximum current is close to that obtained by a much 
more demanding microscopic calculation based on BdG 
equations ( 8, 9]) 

It must be said that in Ref. [8, 9] the vortex struc- 
ture is imposed on the condensate order parameter A. 
Therefore near the core the current goes to zero because 
A vanishes and the superfluid density vanishes with it, 
since ps oc A^. However in the core region, where the 
superfluid velocity exceeds the critical value for the cre- 
ation of single particle excitations, ps < p, so that the 
total density is non zero even at the core where ps = 0- 

On the contrary, in our approach which does not take 
into account single particle excitations, all the fluid is 
superfluid, and therefore the vanishing at the core of the 
superfluid current due to vorticity implies that the total 
density is zero there. 

Results similar to those reported above were obtained 
from the simulations of the lower density system, the 
main difference being the shallower vortex cores (which 
scale as kp^). The calculated superfluid current density 
around the core of an isolated vortex, plotted as a func- 
tion of kpr, is indistinguishable from the one shown in 
Fig. ©. 

We show in Fig. ^ the calculated drag force on 
the object, as a function of the velocity V. Each value 
is obtained as a time-average j2S] of curves like the 
one shown in Fig. © (the average is taken over a 
time interval where the drag force has already reached 
a plateau). From these results, a value of the critical 
velocity Vc ^ 0.4 ci is obtained, in agreement with the 
simple estimate ([S]). 
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FIG. 7: Average drag force for the high and low density sys- 
tems, plotted as a function of V. 

We also find that the drag force first increases linearly 
with the velocity ("Stoke's law", usually associated with 
laminar drag) and then turns to a quadratic behavior 
("Newton's law" , usually associated with turbulent flow), 
since at supersonic velocities there is also a contribution 
to the drag associated with sound radiation. A similar 
behavior is observed in the case of an impenetrable cylin- 
der moving inside a dilute BEC ^23\. 

The phenomenology of the dissipative motion of an 
object displaced through a UFG, as it appears from our 
calculations, is qualitatively similar in many aspects, in 
spite of the different non-linear interactions, to the be- 
havior observed in the case of an object moving in a BEC 
[20I [25| : the occurrence of vortex emissions in pairs and 
the associated density patterns are similar in the two 
systems, and also the behavior of the drag as a function 
of the object velocity. There are important differences 



though: the emitted vortices in the UFG are doubly 
quantized, as expected from fermion pairing; the pre- 
dicted critical velocity is different; the vortex core struc- 
ture is also different, scaling in the present case as the 
inverse Fermi momentum. 



V. CONCLUSIONS 

In conclusion, we have numerically studied the mo- 
tion of an object in the ultracold unitary Fermi gas. We 
described the system by using an extended density func- 
tional approach, which has been used recently to success- 
fully describe a number of static and dynamical proper- 
ties of cold Fermi gases. 

We find that quantized vortices are spontaneously gen- 
erated in pairs during the time evolution at supercritical 
velocities. Moreover, the profile of the current density as 
a function of the distance from the core is quantitatively 
close to the one found by the much more demanding so- 
lution of the BdG solutions. We explain this agreement 
by observing that at distances larger than the one corre- 
sponding to the maximum current density, in both treat- 
ments the superfluid density coincides with the total den- 
sity, since the speed of the fluid is below its critical value 
for the creation of single particle excitations. At shorter 
distances, in the core region, the current density is again 
similar in the two treatments: in both it vanishes imply- 
ing that the superfluid density is zero at the center of 
the vortex. The main difference is that while in our case 
the superfluid density coincides with the total density, in 
the BdG treatment the density may have a contribution 
from a " normal" component related to single particle ex- 
citations and this component provides a nonzero density 
at the center. 

We acknowledge Fondazione CARIPARO (project of 
excellence "Macroscopic Quantum Properties of Ultra- 
cold Atoms under Optical Gonfinement" ) and University 
of Padova (progetto di ateneo "Quantum Information 
with Ultracold Atoms in Optical Lattices") for partial 
support. 
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